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Abstract 

We analyze the asymptotic stability of a nonlinear system of two differential equa- 
tions with delay describing the dynamics of blood cell production. This process takes 
place in the bone marrow where stem cells differentiate throughout divisions in blood 
cells. Taking into account an explicit role of the total population of hematopoietic stem 
cells on the introduction of cells in cycle, we are lead to study a characteristic equation 
with delay-dependent coefficients. We determine a necessary and sufficient condition for 
the global stability of the first steady state of our model, which describes the population's 
dying out, and we obtain the existence of a Hopf bifurcation for the only nontrivial pos- 
itive steady state, leading to the existence of periodic solutions. These latter are related 
to dynamical diseases affecting blood cells known for their cyclic nature. 

Keywords: asymptotic stability, delay differential equations, characteristic equation, delay- 
dependent coefficients, Hopf bifurcation, blood cell model, stem cells. 

1 Introduction 

Blood cell production process is based upon the differentiation of so-called hematopoietic 
stem cells, located in the bone marrow. These undifferentiated and unobservable cells have 
unique capacities of differentiation (the ability to produce cells committed to one of the three 
blood cell types: red blood cells, white cells or platelets) and self-renewal (the ability to 
produce cells with the same properties). 

Mathematical modelling of hematopoietic stem cells dynamics has been introduced at 
the end of the seventies by Mackey [5T]. He proposed a system of two differential equations 
with delay where the time delay describes the cell cycle duration. In this model, hematopoi- 
etic stem cells are separated in proliferating and nonproliferating cells, these latter being 
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introduced in the proliferating phase with a nonlinear rate depending only upon the non- 
proliferating cell population. The resulting system of delay differential equations is then 
uncoupled, with the nonproliferating cells equation containing the whole information about 
the dynamics of the hematopoietic stem cell population. The stability analysis of the model 
in [21] highlighted the existence of periodic solutions, through a Hopf bifurcation, describing 
in some cases diseases affecting blood cells, characterized by periodic oscillations [TT?] . 

The model of Mackey [5T] has been studied by many authors, mainly since the beginning 
of the nineties. Mackey and Rey [531 HH [25] numerically studied the behavior of a struc- 
tured model based on the model in [21], stressing the existence of strange behaviors of the 
cell populations (like oscillations, or chaos). Mackey and Rudnicky (2S1 [27] developed the 
description of blood cell dynamics through an age-maturity structured model, stressing the 
influence of hematopoietic stem cells on blood production. Their model has been further 
developed by Dyson et al. [HI UU [T5] , Adimy and Pujo-Menjouet [7J, Adimy and Crauste 
[21 |3J and Adimy et al. [I] . Recently, Adimy et al. [5J E] studied the model proposed in 
[2"T] taking into account that cells in cycle divide according to a density function (usually 
gamma distributions play an important role in cell cycles durations), contrary to what has 
been assumed in the above-cited works, where the division has always been assumed to occur 
at the same time. 

More recently, Pujo-Menjouet and Mackey [3(1] and Pujo-Menjouet et al. [21] gave a 
better insight into the model of Mackey [21] , highlighting the role of each parameter of the 
model on the appearance of oscillations and, more particularly, of periodic solutions, when 
the model is applied to the study of chronic myelogenous leukemia [16) . 

Contrary to the assumption used in all of the above-cited works, we study, in this paper, 
the model introduced by Mackey ;2T] considering that the rate of introduction in the prolifer- 
ating phase, which contains the nonlinearity of this model, depends upon the total population 
of hematopoietic stem cells, and not only upon the nonproliferating cell population. The in- 
troduction in cell cycle is partly known to be a consequence of an activation of hematopoietic 
stem cells due to molecules fixing on them. Hence, the entire population is in contact with 
these molecules and it is reasonable to think that the total number of hematopoietic stem 
cells plays a role in the introduction of nonproliferating cells in the proliferating phase. 

The first consequence is that the model is not uncoupled, and the nonproliferating cell 
population equation does not contain the whole information about the dynamics of blood 
cell production, contrary to the model in [21, 29, 30J. Therefore, we are lead to the study of 
a modified system of delay differential equations (system ©-([I])), where the delay describes 
the cell cycle duration, with a nonlinear part depending on one of the two populations. 

Secondly, while studying the local asymptotic stability of the steady states of our model, 
we have to determine roots of a characteristic equation taking the form of a first degree 
exponential polynomial with delay-dependent coefficients. For such equations, Beretta and 
Kuang [5] developed a very useful and powerful technic, that we will apply to our model. 

Our aim is to show, through the study of the steady states' stability, that our model, 
described in ([HI)-®, exhibits similar properties than the model in [2T] an d that it can be 
used to model blood cells production dynamics with good results, in particularly when one is 
interested in the appearance of periodic solutions in blood cell dynamics models. We want to 
point out that the usually accepted assumption about the introduction rate may be limitative 
and that our model can display interesting dynamics, such as stability switches, that have 
never been noted before. 

The present work is organized as follows. In the next section we present our model, 
stated in equations |3j and ([4]). We then determine the steady states of this model. In 
section [33 we linearize the system ©-Q about a steady state and we deduce the associated 
characteristic equation. In section 01 we establish necessary and sufficient conditions for the 
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global asymptotic stability of the trivial steady state (which describes the extinction of the 
hematopoietic stem cell population). In section [5l we focus on the asymptotic stability of the 
unique nontrivial steady state. By studying the existence of pure imaginary roots of a first 
degree exponential polynomial with delay-dependent coefficients, we obtain the existence of 
a critical value of the time delay for which a Hopf bifurcation occurs at the positive steady 
state, leading to the appearance of periodic solutions. Using numerical illustrations, we show 
how these solutions can be related to periodic hematological diseases in section [SJ and we 
note the existence of a stability switch. We conclude with a discussion. 

2 A Nonlinear Model of Blood Cell Production 

Let consider a population of hematopoietic stem cells, located in the bone marrow. These 
cells actually perform a succession of cell cycles, in order to differentiate in blood cells (white 
cells, red blood cells and platelets). According to early works, by Burns and Tannock [TT] 
for example, we assume that cells in cycle are divided in two groups: proliferating and 
nonproliferating cells. The respective proliferating and nonproliferating cell populations are 
denoted by P and N. 

All hematopoietic stem cells die with constant rates, namely 7 > for proliferating cells 
and 6 > for nonproliferating cells. These latter arc introduced in the proliferating phase, in 
order to mature and divide, with a rate (3. At the end of the proliferating phase, cells divide 
in two daughter cells which immediately enter the nonproliferating phase. 

Then the populations P and N satisfy the following evolution equations (see Mackey [ST] 
or Pujo-Menjouet and Mackey [3"0]). 



In each of the above equations, r denotes the average duration of the proliferating phase. The 
term e _7T then describes the survival rate of proliferating cells. The last terms in the right 
hand side of equations |JJ and account for cells that have performed a whole cell cycle and 
leave (enter, respectively) the proliferating phase (the nonproliferating phase, respectively). 
These cells are in fact nonproliferating cells introduced in the proliferating phase a time r 
earlier. The factor 2 in equation (0) represents the division of each proliferating hematopoietic 
stem cell in two daughter cells. 

We assume that the rate of introduction (3 depends upon the total population of hematopoi- 
etic stem cells, that we denote by S. With our notations, S = P+N. This assumption stresses 
the fact that the nature of the trigger signal for introduction in the proliferating phase is the 
result of an action on the entire cell population. For example, it can be caused by molecules 
entering the bone marrow and fixing on hematopoietic stem cells, activating or inhibiting 
their proliferating capacity. This occurs in particularly for the production of red blood cells. 
Their regulation is mainly mediated by an hormone (a growth factor, in fact) called erythro- 
poietin, produced by the kidneys under a stimulation by circulating blood cells (see Belair 
et al. [5], Mahaffy et al. [25] )• 

Hence we assume that 



The function (3 is supposed to be continuous and positive on [0, +00), and strictly decreasing. 
This latter assumption describes the fact that the less hematopoietic stem cells in the bone 




jP{t) + (3N{t) - e~ 1T f3N{t - t), 



SN(t) - (3N(t) + 2e^ T pN(t - r). 



(1) 



(2) 
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marrow, the more cells introduced in the proliferative compartment |21l 130] . Furthermore, 
we assume that 

lim (3{S) = 0. 

Adding equations ([T]) and we can then deduce an equation satisfied by the total 
population of hematopoietic stem cells S(t). We assume, for the sake of simplicity, that 
proliferating and nonprolifcrating cells die with the same rate, that is S = 7. Then the 
populations N and S satisfy the following nonlinear system with time delay r, corresponding 
to the cell cycle duration, 

^(t) - -5S(t) + e-^(3(S(t-T))N(t-r), (3) 

^(t) = ~5N(t)-P(S(t))N(t) + 2e- s Tp(S(t-T))N(t-T). (4) 

From Hale and Verduyn lunel 18J, for each continuous initial condition, system ©-((I]) 
has a unique continuous solution (S(t), N(i)), well-defined for t > 0. 

Lemma 2.1. For all nonnegative initial condition, the unique solution {S(t), N(t)) of ©l-Q) 
is nonnegative. 

Proof. First assume that there exists £ > such that N(£) = and N(t) > for t < £. 
Then, from Q and since /? is a positive function, 

^(0=2e-^(5(^-T))iV(C-T)>0. 

Consequently, iV(t) > for t > 0. 

If there exists C > such that S(() = and S(t) > for t < £, then the same reasoning, 
using ((3]), leads to 

^(C)=e- 5 ^(5(C-r))iV(C-r)>0, 

and we deduce that S(t) > for t > 0. □ 

Remark 1. The positivity of S and N, solutions of system does not a priori implies 

that P = S — N is nonnegative. 

Using a classical variation of constant formula, the solutions P(t) of (QJ) are given, for 
t > 0, by 

P(t) = e- St P(0) + e' st [ e se l3{S(0))N(9) - e &{9 ~ r) /3{S(6 - t))N(9 - r)d9. 
Jo 

Setting the change of variable a — 9 — t, we obtain 

r r° 

P{t) = e- st P(0) - J e se /3(S(9))N{9)d9 
Consequently, P(t) > for t > if 

r° 

P(0) > J e se f3(S(9))N(9)d9, 

that is if 

r° 

S(0) > N(0) + / e S9 /3(S(9))N(9)d9. 



t 

e- 5t I e se f3(S(9))N{9)d9. 

t-T 
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This condition is biologically relevant since f® T e se /3(S(9))N(6)d6 represents the population 
of cells that have been introduced in the proliferating phase at time 9 S [— r, 0] and that 
have survived at time t = 0. Hence, from a biological point of view, the population in the 
proliferating phase at time t — should be larger than this quantity. 

In the stability analysis of system I©-©, the existence of stationary solutions, called 
steady states, is relevant since these particular solutions are potential limits of system ([3])- 
®. ' __ 

A steady state of system ©-© is a solution (S, N) satisfying 

dS _ d~N _ 
dt ~ dt ~ 

Let (S,N) be a steady state of ©-©. Then 

SS = e- Sr (3(S)N, (5) 
(5 + (3(S))N = 2e~ 5T {3(S)N. (6) 

We immediately notice that (0, 0) is a steady state of system ©-© , that we will denote, 
in the following, by E°. This steady state always exists. It describes the extinction of the 
hematopoietic stem cell population. 

Assume that (S,N) is a steady state of ©-© with S,N ^ 0. Then, from © and @, 

(2e- 5T - 1)/3(S) = 6 and N = 



e- ST f3(S) ' 

A necessary condition to obtain a nontrivial steady state is then 2e~ dT — 1 > 0, that is 

r< ln(2)_ 
o 

Under this condition, since the function (3 is decreasing, positive, and tends to zero at infinity, 
there exists S > satisfying 

(2e- 5T - l)p(S) = 5, (7) 

if and only if 

(2e- ST - 1)0(0) > S. (8) 
One can easily check that condition ([S]) is equivalent to 

6< m and 0<.<T:=im(A). (9) 

In this case, S > solution of ((7|) is unique, and 

e -8 T 

One can check that N < S. These results are summed up in the next proposition. 

Proposition 2.1. // inequality holds true, the system f3p-|7p has exactly two steady 
states: E° = (0,0) and E* — (S*,N*), where S* > is the unique solution of equation |7p 
and N* = (2e- 5r - l)e 5r S*. 
If 

[2e- 5T - 1)13(0) < S, 
then system has only one steady state, namely E° = (0,0). 

In the next section, we linearize the system ©-((H) about one of its steady states in order 
to analyze its local asymptotic stability. 
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3 Linearization and Characteristic Equation 

We are interested in the asymptotic stability of the steady states of system To that 

aim, we linearize system about one of its steady state and we determine the associated 

characteristic equation. We assume that is continuously differentiable on [0, +00). 

Let (S, N) be a steady state of system ©-Q- From Proposition (S,N) is either E° 
orE*. 

The linearization of system ©-(H]) about (£, N) leads to the following system, 

7 Q 

— (t) = -SS(t) + e- ST f3(S)N(t-T)+e- ST N(3'(S)S(t-T), (10) 
dt 

+2e- gT [P{S)N(t - T )+ NP'(S)S(t - t)] , (11) 

where we have used the notations S(t) and N(t) instead of S(t) — S and N(t) — N for the 
sake of simplicity. 

The system ([TU)) - ((TT1) can be written 

( m ) + A( S{t - T) ) 



where 



-5 \ / a P(S) 

1 -a -(S + p(S)) I - \ 2a 2(3{S) 



and 

a = a(N,S) :=N/3'(S). (12) 

The characteristic equation of system (flQ|) -([TT j) associated with the steady state (S, N) is 
defined by 

det(A -A\- e- XT A 2 ) = 0. 
After calculations, this equation reduces to 

(A + 5) [A + 5 + f3(S) - (2{3(S) + a(N, S))e- Sr e~ Ar ] = 0. (13) 



We recall that the steady state (S, N) is locally asymptotically stable when all roots of (TH 
have negative real parts and the stability can only be lost if eigenvalues cross the imaginary 
axis, that is if pure imaginary roots appear. 

One can notice that A = — 5 < is always an eigenvalue of (fTS"]) . Therefore, we only focus 
on the equation 

A + 8 + /3(S) - (2/3(S) + a(N, S))e- 5T e - XT = 0. (14) 

We first analyze, in the next section, the stability of the trivial steady state E°. We 
establish necessary and sufficient conditions for the population's dying out. Then, in section 
03 we concentrate on the behavior of the positive steady state E* . 
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4 Global Asymptotic Stability of the Trivial Steady State: 
Cell's Dying Out 

We concentrate, in this section, on the stability of the steady state E° = (0,0). From (fT2"|) . 
a(0, 0) = 0, so, for S = N = 0, the characteristic equation (TT4"|) becomes 

A + 5 + /3(0) - 2(3{Q)e- ST e - XT = 0. (15) 

It is straightforward to see that equation (|15p has a unique real eigenvalue, say Ao, and all 
other eigenvalues A ^ A of (fTS"]) satisfy Re(A) < Ao- 

Let consider the mapping A ^ A + <5 + 0(0) — 20(O)e~ Ar e _AT as a function of real A. Then 
it is an increasing function from — oo to +oo, yielding the existence and uniqueness of Ao- 

Assume that A = /i + ioj ^ Ao satisfies (fTS")) . Then, considering the real part of ([15]) , we 

get 

At - Ao = 2/3(0)e- Sr [e-^ cos(wr) - e - XaT ] . 

By contradiction, we assume that /i > Ao- Then e~' IT cos(cjt) — e~ A ° r < and we obtain a 
contradiction. So /x < Ao- Now if \x = Ao, then the previous equality implies that 

cos(wr) = 1, for r > 0. 

It follows that sin(ojr) = and, considering the imaginary part of (|15p with A = /i + iuj, 
given by 

u + 2[3{Q)e- ST e-> iT sin(wT) = 0, 

we obtain u> = and A = Ao, which gives a contradiction. Therefore, fi < Xq. 
The real root Ao is negative if 

[2e- &T - l)/3(0) < 6, 

and all eigenvalues of (fTS"]) have negative real parts in this case. When condition |8]) holds, 
Ao is positive. We can then conclude in the next proposition to the stability of E°. 

Proposition 4.1. The trivial steady state E° = (0,0) of system f3f)-|^P is locally asymptot- 
ically stable when 

{2e~ ST - 1)0(0) < 5, (16) 

and unstable when 

(2e- ST - 1)0(0) > S. 

Remark 2. When 

(2e~ ST - 1)0(0) = 6, 

then the unique real root of \15\) is Ao = 0, and all other eigenvalues have negative real parts. 
One can easily check that Ao = is a simple root of \15\) , since the first derivative of the 
mapping AhA + <5 + (3(0) - 2/3(0)e _5T 'e _Ar at A = is 

1 + 2f3{0)re- ST > 0. 

Then the linear system is stable, but we cannot conclude to the asymptotic stability of the 
trivial steady state E — (0, 0) of system (0)-0) without further analysis. This is done in 
Proposition \4-'A 

When condition (fT6|) holds true, E° is in fact the only steady state of system ((3])-((4|) (see 
Proposition ^. ip . In this case, we can show that E° is globally asymptotically stable. 
We first show the following result. 
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Lemma 4.1. Let (S(t), N(t)) be a solution of If 'lim t _> +00 N(t) = 0, then lim t ^ +00 S(t) 
0. 

Proof. Using ([3]), a classical variation of constant formula gives us, for t > 0, 

S(t) = e- St S(0) + e- st f e 5{e ~ T) /3(S(6 - t))N(8 - r)d6. 
Jo 

Setting ci — — t, this expression becomes 



S{t) = e~ dt 5(0) + 



-61 



5a f3(S(a))N(a)da. 



Let e > be fixed. Since N is assumed to tend to zero when t tends to oo, there exists T > 
such that 

N{t)<e-^-;, fort>T. (17) 



2/3(0) ' 



Then, for t > T + r, 



S(t) = 



,-<5t 



5(0)+ / e 5,y f3(S(a))N(a)da 



-Si 



J f3(S(a))N(a)da. 



Using (TT71) an d the fact that /3(0) is a bound of /3, we obtain, for t > T + r, 



S(t) < e 



< e 



< e 



-5t 



-<5t 



5(0) 
5(0) 
5(0) 



J f3(S(a))N(a)da 



e Sa p{S{a))N(a)da 



V(--* 



e B °P{S{(T))N((T)d<T 



Let t > be such that 



-5* 



5(0)+ / e Sa p{S{a))N{o-)do- 



<2' 



for i > t 



Then, for t > max{i, T + r}, we obtain 

5(t) < e. 

Thus, 5(f) tends to zero as t tends to +oo, and the proof is complete. 

We recall a very useful lemma, proved by Barbalat (see Gopalsamy [T7] 



□ 



Lemma 4.2. Let f : [a, +oo) — > R ; a G R, fee a differ entiable function. If lim t ^ +OQ /(i) 
exists and f'(t) is uniformly continuous on (a,+oo), i/ien 

lim f'(t) = 0. 



We then prove the following result, dealing with the global asymptotic stability of E°. 
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Proposition 4.2. Assume that 

(2e~ ST - 1)0(0) < 5. (18) 

Then all solutions (S(t), N(t)) of system 0]-^]) converge to the trivial solution (0, 0). Hence 
E° is globally asymptotically stable and the cell populations dye out. 

Proof. Let (S(t),N(t)) be a solution of We define, for t > 0, 



Y(t) = N(t) + 2e- Sr I 0(S(6))N(9)dd 



t-T 

Using (3]), we can check that 

Y'(t) = N(t) [{2e~ ST - l)0(S(t)) - S] . (19) 
From condition (fT5|) and the fact that is decreasing, it follows that 

Y'(t) < for t > 0. 

Thus, Y is decreasing. Since Y is a nonnegative function, we deduce that there exists y > 
such that 

lim Y(t) = y. 

t — >+oo 

In particularly, Y is bounded, and consequently N is also bounded. 

We then deduce, with (HI, that N' is bounded and, using a similar technic than the one 
used in the proof of Lemma |4~T1 that S is bounded. Consequently, with ©, we obtain that 
S' is bounded. 

From (fT9|) . since N, S, N' and S' are bounded, Y' is uniformly continuous. 
Since lim t ^ +00 Y(t) exists and Y' is uniformly continuous on (0, +oo), Lemma implies 
that 

lim Y'(t) = 0. 

t — > + OG 



Consequently, from IT9|) . we obtain either 



lim N(t) = or lim {2e~ aT - l)(3(S(t)) = 5. 

t — >+oo t — >-|-oo 



First, assume that (TIT))) holds true, that is 

(2e- 5r - 1)0(0) < 5. 

If 2e~ bT — 1 > 0, then, since (3 is a decreasing function satisfying (TH)|) , we deduce that 
(26" 157 " - l)0(S(t)) < (2e- Sr - 1)13(0) < S. If 2e~ Sr - 1 < 0, then (2e- 5r - 1)0 (S(t)) <0<6. 
Consequently, if it exists, lim t ^ +00 (2e _l5T — l)f3(S(t)) cannot be equal to S and it follows 
that 

lim N(t) = 0. 

t — H-co 

From Lemma |4~T1 we deduce that lim t ^ +00 S(t) = 0, and the conclusion follows. 
Second, assume that 

(2e- 5T - 1)0(0) = S. 

Then, lim t ^ +00 (2e~ 5r - l)(3(S(t)) = S is equivalent to lim t ^ +00 /3(S(t)) = 0(0). Since is 
positive and decreasing, this is equivalent to lim t ^ +00 S(t) = 0. It follows that either 

lim N(t) =0 or lim S(t) = 0. 
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If lim t ^ +oc N(t) = 0, we conclude similarly to the previous case with Lemma 14.11 So we 
assume that lim f ^ +00 S(t) = 0. 
From ©, we deduce that 

lim (3(S(t - r))N{t - t) = 0. 

t — >+oc 

Consequently, either 

lim P(S(t - r)) = or lim N(t - r) = 0. 

t— » + oo t—>-\-oo 

Since lim t ^ +oc S(t) = 0, then lim f ^ +oc f3(S(t-r)) = /3(0) > 0. Hence, lim t ^ +00 N(t-r) = 0, 
and it follows that lim t ^ +oc N(t) — 0. 

This concludes the proof. □ 

Remark 3. Let C denote the set of continuous functions mapping [— r, 0] into R + . One can 
check that the function V , defined for (ip, ip) £ C 2 fry 

^) - ^(0) + 2e~ Sr J (3{ V {6))^{0)dO, 

satisfies 

V{tp, 1>) = V(0) [(2e~^ - l)/%(0)) - <J] . 
Hence, V is a Lyapunov functional (see Hale and Verduyn Lunel JTE/) on the set 

G= {(^) eC 2 ; ^(0) [(2e~ 4T - l)/3(^(0)) - <S] < 0} . 

W^/i assumption tlfy) , G = C 2 . In the proof of Proposition \JJ^ we did not directly use the 
properties of Lyapunov functionals, but the function Y is defined by 

Y(t)=V(S t ,N t ), fort>0, 

where St (respectively, N t ) is defined by St (6) = S(t + 9) (respectively, N t {9) — N(t + 6)), 
0€[-t,O]. 

Through Propositions 14.11 and 14. 2\ we obtained necessary and sufficient conditions for 
the global asymptotic stability of E°. Therefore, in the next section, we concentrate on the 
behavior of E*, the unique positive steady state of ©-([1]). 

5 Local Asymptotic Stability of the Positive Steady State 

We now turn our considerations on the stability of the unique nontrivial steady state of system 
©-(H), namely E* = (S*,N*), where, from Proposition ^. 11 S* is the unique solution of Q, 
and TV* = (2e~ ST - l)e Sr S*. 

In order to ensure the existence of this steady state, we assume that condition l[5]). or 
equivalently condition holds true. That is 

6< m ,„d 0< T < T :.Il„( J »_)^ 

In particularly, 2e~ bT — 1 > 0. 

In this case, Proposition 14. II indicates that the unique other steady state E° is unstable. 
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From their definitions in Proposition 12. 1[ the steady states S* and N* depend on the 
time delay r. In fact, 

S* = S*(T)=f3- 1 (^J r ^) and N* = iV*(r) = ^zl S *(r), 

where /3" 1 : (0, [3(0)] — > [0, +oo) is a decreasing function. 

Using these expressions, we can stress that S* and N* are positive decreasing continuous 
functions of r £ [0,r), continuously differentiable, such that S*(0) = N*(0) — /3 _1 (<5), and 
]im T ^(S*(T),N*(r)) = (0,0) = E°. _ _ 

The characteristic equation (fTi)) . with S = S* and AT = TV*, is then given by 

A + A(t) — B(t)c~ Xt = 0, (20) 

with 

A(t) :=S + (3(S*(t)) and B(t) := [2f3(S*(r)) + N*(T)f3'(S*(T))]e- ST . 

Notice that A(t) > for all r S [0,r). Moreover, from (|7|), we obtain 

B(t) =A(T) + (2e- ST -l)S*(T)/3'{S*{T)), forre[0,r). (21) 

In particular, B(t) < A{r) for r e [0,r). 
Taking r = in PU| . we obtain 

A + A(0) - 5(0) = 0, 

that is 

Since /? is decreasing, we deduce that the only eigenvalue of PD]) is then negative. The 
following lemma follows. 

Lemma 5.1. When 5 < /3(0) and r = 0, the nontrivial steady-state E* of system f3Jj-|^P is 
locally asymptotically stable, and the system undergoes a transcritical bifurcation. 

When r increases and remains in the interval [0,t), the stability of the steady state can 
only be lost if purely imaginary roots appear. Therefore, we investigate the existence of 
purely imaginary roots of (f20|) . 

Let A = iui, lo £ K, be a pure imaginary eigenvalue of (f^U|) . Separating real and imaginary 
parts, we obtain 

A(t) - B{t) cos(wr) = 0, (22) 
lo + S(r)sin(wr) = 0. (23) 

One can notice, firstly, that if u> is a solution of f2"2"|) - ((2"3"|) then — uj also satisfies this system. 
Secondly, to — is not a solution of f2"2"|) - ((2"3"|) . Otherwise, we would obtain A(t) = B(t) for 
some r £ [0,r), which contradicts £?(t) < v4(t) for r £ [0,r). Therefore ui — cannot be 
a solution of (l22 |) -([23 |) . Thus, in the following, we will only look for positive solutions uj of 

(EH)-®. 

From (|22[) . a necessary condition for equation (|20p to have purely imaginary roots is that 

A(r) < |B(r)|. 
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Since B(t) < A(t), this implies in particularly that B(t) must be negative. Moreover, from 
([2"Tj) . the above condition is equivalent to 

2A(r) + (2e- 5r - 1)S*(t)/3'(S*(t)) < 0. 

Using the definitions of A(t), N*(t), S*(t) and equality this inequality becomes the 
following condition on r, 



2e- 5r -l y ' \2e- ST -lj' \ \2 e - s 

Lemma 5.2. Let x '■ [0, +oo) — > (— oo,0] &e defined, for y >0, by 

xiv) = vP'iv)- 

Assume that 

(Hi) x * s decreasing on the interval [0,/3 _1 ((5)] . 
(H 2 ) X (0-HS))<-M. 

Then there exists a unique r* G (0, r) such that condition \2J$ is satisfied if and only if 
re[0,-r*). 

Proof. Consider the negative functions /i(r) and / 2 (t), defined for r G [0, r] by 



/xW=x(r 1 (^-r')) and / a (r) = - 4 *' 



2e" 5T - I// ' JAy ' {2e- ST - l) 2 ' 

Then 

{r G [0,t) ; condition ([24]) is satisfied} = {r G [0,r) ;/i(r) < / 2 (r)} . 
The function f 2 satisfies 

/ 2 (0) = -4«5 and f 2 (r) = -2 S + f ( ° } ,3(0), 



/aW = ^—57 775 < °- 



and, for r G [0, r], 

4(5 2 e- 5r (2e- 5r + 1 
(2e- 5T - 1) 

Hence / 2 is decreasing from —AS to — 2(6 + f3(0))[3(0)/S. 
Forre[0,T), 

Since is decreasing on (0,/3(0)] and, from (Hi), x is decreasing on [0,/3 _1 (5)], we deduce 
that /i is increasing. Moreover, 

/i(0)=x(/3- 1 («5)) and /i(r) = 0. 

From (H 2 ), 

/i(0)</ 2 (0). 

Since /i(r) > / 2 (t), there exists r* G (0,r), which is unique since f\ is increasing and / 2 
decreasing, such that 

{re[0,r) ;/i(r)</ 2 (r)}=[0,r*). 
This concludes the proof. □ 
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Remark 4. Assumption (Hi) is not necessary for the existence of t* , it just implies the 
uniqueness. This latter can be achieved with weaker conditions, as we will check on an 
example in section^ 

We assume, in the following, that (Hi) and (H2) are fulfilled and r £ [0, r*). 
System (j2"2"|) - (|2"3"l) is equivalent to 

1 \ A(t) . ui . s 

cos(cjr) = bJt) 1 sm( ] = ~W)- ( } 

Note that for r e [0, r*), B(t) < 0. 

Therefore, adding the squares of both sides of (f25|) . purely imaginary eigenvalues iui of 
(|2T))) . with a; > 0, must satisfy 

w = VS 2 (r)-4 2 (r). (26) 

So, in the following, we will think of ui as w(t). 

Substituting this expression for in (|25p . we obtain 



B(t)' 

(27) 



:os (rv/B2(r)-A2(r)) = f 7) 
sin(V ^(r)-^(r)) - V^W^(7 



From the above reasoning, values of r 6 [0,r*) solutions of system (|27f generate positive 
uj(t), given by (|2l)|). and hence yield imaginary eigenvalues of (|2H)) . Consequently, we look 
for positive solutions r of l|27p in the interval (0, r*). 
Positive solutions t 6 (0,t*) of (|2T]) satisfy 

t^B 2 (t) - A 2 (t) = arccos ( 4t4 ) + 2/br > ^N , 

\ S ( T )/ 

where No denotes the set of all nonnegative integers. We set 

arccos ( 4£4 ) + 2kn 
r k (r) = ) ( U= - k E N ,r £ [0,r*). 



Values of t for which w(t) = J B 2 (t) — A 2 (t) is a solution of are roots of the functions 

z k (T) = T-T k ( T ), ken ,r S [0,T*). (28) 

The roots of can be found using popular software like Maple, but are hard to determine 
with analytical tools [9]. The following lemma states some properties of the Z k functions. 

Lemma 5.3. For k E No, 

Z k (0) < and lim Z k (r) = -00. 

T — >T * 

Therefore, provided that no root of Z k is a local extremum, the number of positive roots of 
Z k , k £ No, on the interval [0, r*) is even. 

Moreover, if Z k has no root on the interval [0, t*), then Zj, with j > k, does not have 
positive roots. 
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Proof. Notice first that Tfe(O) > and, secondly, that 

f A(t)\ 

lim arccos ( I = 7r and lim \J B 2 (r) — A 2 {t) = 0, 

so lim T ^ r » Tfc(r) = +00. Then the first statement holds. 

To prove the second statement, one can notice that Tk+i{r) > Tfc(r), for r G [0, r*) and 
k e N , so 

Z k+1 {r) < Z k {r), re[0,r'),ieN„. 
Using the fact that Z k (0) < 0, we conclude. This ends the proof. □ 



Remark 5. The second statement in Lemma 15.31 implies, in particularly, that, if Zq has 
no positive root, then \25)) has no positive solution, and equation \20\) does not have pure 
imaginary roots. 

In the following proposition we establish some properties of pure imaginary roots of 
equation ((2"U|) . using a method described in [30]. 

Proposition 5.1. Let ±iaj(r c ), with lo(t c ) > 0, be a pair of pure imaginary roots of equation 
when t — t c . Then Hlo(t c ) are simple roots of \20}) such that 



sign < 



dRe{\) 



dr 



signl — B — B 2 B't c + B(A + A' + AA't c ) - B'A } , (29) 



where A = A{t c ), B = B(t c ), A' = A'(t c ) and B' = B'(t c ). 
Proof. We set 

A(A, t) = A + A(t) — B{r)e~ XT . 

Let A(r) be a family of roots of (j2"D|) . so A(A(t), r) = 0, such that A(t c ) is a pure imaginary 
root of ([2"U)l . given by A(r c ) = iuj(T c ). Then, 

^(r)A A (A,r) + A r (A,r) =0, (30) 

where 

Aa(A,t) :=^(A,r) = l + B(r)re- AT , 

and 

A T (A,r) := ^(A,r) = A'(r) + [B(r)X - B'(r)]e- X \ 

Assume, by contradiction, that A(r c ) = ioj(T c ) is not a simple root of (|2"0|) . Then, from 
CT. A T (iw(r c ),r c ) = A'(r c ) + [iS(r c )w(r c ) - S'(r c )]e-^( T <=)' rc = 0. Separating real and 
imaginary parts in this equality we deduce 

B'{t c )cos{lj(t c )t c ) - B(T c )u(T c )sm(uj(T c )T c ) = A'(t c ), 
B(t c )uj(t c )cos(uj(t c )t c ) + B'(T c )sm(uj(T c )T c ) = 0. 

We recall that B(t c ) is necessarily strictly negative. Using (|2"5|) , the above system is equivalent 
to 

a _ A(r c )B'(r c ) - S(r c )A'(r c ) 



w(t c ) 



S'(r„)- 
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Since w(r c ) > and satisfies, from uj(t c ) 2 = B(t c ) 2 — A(t c ) 2 , we obtain 

A(t c )B'(t c ) - B{t c )A'{t c ) 



B(t c ) — A(t c ) = 

A(t c )B(t c ) = -B'(t c ). 



B(t c ) 



Substituting the second equation in the first one, this yields B(t c ) 2 — A{t c ) 2 = — A'(t c ) — 
A(t c ) 2 , so 

B(t c ) 2 = -A'(t c ). 

Since B(t c ) 2 > and A'(t c ) > 0, we obtain a contradiction. Hence A(t c ) — ioj(r c ) is a simple 
root of (f2U|> . 

In the following, we do not mention the dependence of the coefficients A and B (and their 
derivatives) with respect to r. 



Now, from (1501) , we obtain 



Since A(A, r) = 0, we deduce 



e XT + Bt 



B' - BX - A'e XT ' 



B 



Therefore, 



For t = t>, we obtain 



Then, 



Re 



dX 



dT 

Noticing that 



X + A 

1 _ B + Bt(X + A) 
~ (B' - BX)(X + A) - A'B ' 

B + BT c (ioj(T c ) + A) 
(B' - iBuo{T c )){iuo{T c ) + A) - A'B ' 

B(l+AT c )+iBuj{T c )T c 

B'A - AB' + Boj 2 (t c ) + i(B' - AB)uj{t c ) ' 

[B 2 (l + At c ) + Bt c {B' - AB)]u{t c ) 2 + B(l + At c ){B'A - A'B) 
[B'A - A'B + Boj(t c ) 2 } + [B' - AB} 2 oj(t c ) 2 ' 



. f dRe(X) 1 



hi;; It { Re | ^ 



we get, 

sign | | = sign^[B 2 (l + AT c )+BT c (B' -AB)]lu(t c ) 2 

+B(1 + At c )(B'A - A'B) 

Since zcj(t c ) is a purely imaginary root of (|20|l . then, from 

W (r c ) 2 = B 2 - A 2 . 



(31) 
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Substituting this expression in pip , we obtain, after simplifications 

sign | | = sign js [B 3 + B 2 B't c - B(A 2 + A' + AA't c ) + B'A] J . 

As we already noticed, if equation ([20]) has pure imaginary roots then necessarily B < 0. We 
then deduce (|29|) and the proof is complete. □ 



Using this last proposition and the previous results about the existence of purely imagi- 
nary roots of (|13p , we can state and prove the following theorem, dealing with the asymptotic 
stability of E*. 

Theorem 5.1. Assume that (0J) holds true and (Hi) and (H%) are fulfilled. 

(i) If Zq (defined in i28\) ) has no root on the interval [0,r*), r* defined in Lemma \ 5.2\ 
then the positive steady state E* — (S*,N*) of is locally asymptotically stable 

for t G [0,r). 

(ii) If Zq has at least one positive root t c G (0,t*) then E* is locally asymptotically stable 
for t G [0, r c ) and a Hopf bifurcation occurs at E* for r = r c if 

B(A 2 +A' + AA't c ) - B 3 - B 2 B't c - B'A ^ 0, (32) 

where A = A(t c ), B = B(r c ), A' = A'(t c ) and B' = B'(t c ). 

Proof. First, from Lemma HTT| we know that E* is locally asymptotically stable when r = 0. 

If Zq has no positive root on the interval (0, t*), then the characteristic equation (fT3"|) has 
no pure imaginary root (see Remark [5] and Lemma I5.3|) . Consequently, the stability of E* 
cannot be lost when r increases. We obtain the statement in (i). 

Now, if Zq has at least one positive root, say r c G (0,t*), then equation Q13p has a 
pair of simple conjugate pure imaginary roots ±io;(t c ) for r = r c . From (|3"2"|) together with 
Proposition we have either 



dRe(X) 



dr 



dRe(A) 

> or 



dr 



< 0. 



By contradiction, we assume that there exists a branch of characteristic roots A(r) such that 
A(t c ) = ioj c and 

dRe(X(r)) ^ o 
dr 

for t < t c , t close to r c . Then there exists a characteristic root A(r) such that Re(A(r)) > 
and t < t c . Since E* is locally asymptotically stable when r = 0, applying Rouche's Theorem 
[12] , we obtain that all characteristic roots of (fT3f have negative real parts when r G [0, r c ), 
and we obtain a contradiction. Thus. 



dRe(A) 



> 0. 



dr 

In this case, a Hopf bifurcation occurs at E* when r = t c . □ 

The result stated in (ii) leads, through the Hopf bifurcation, to the existence of periodic 
solutions for system ©-((l]). 

In the next section, we apply the above-mentioned results of stability to a particular 
introduction rate (3 and we present some numerical illustrations. 



16 



F. Crauste 



Stability of a blood cell production model 



6 Example and Numerical Simulations 

We develop, in this section, numerical illustrations of the above mentioned results (mainly 
the ones stated in Theorem 15. ip . 

Let define (see [2H [22 EH IH] ) the introduction rate by 

^ s )=^ ¥rT ^, /M>0, n>l. 

The parameter /3o represents the maximal rate of introduction in the proliferating phase, 9 
is the value for which (3 attains half of its maximum value, and n is the sensitivity of the 
rate of reintroduction. The coefficient n describes the reaction of /3 due to external stimuli, 
the action of a growth factor for example (some growth factors are known to trigger the 
introduction of nonproliferating cells in the proliferating phase [51 \2~E\). 

Then, from ©, the unique positive steady state of ©-([3]) exists if and only if 

5</3o and < r < r = \ In ' 



From ©-([5]), it is defined by 

^ = g( (2e^-l)/ 3 Q_ i y / " ^ N , =e 2e-^l((2e^ 



Note that the function is defined by P' 1 (x) = 9((3 /x - l) 1 /™ for x G (0,/3 ]. 
After computations, we can state that condition (|24p is equivalent to 



i /2AM 

t < t := - In — — — 

5 \ n{0 o + S) 

provided that 

4/3 



Noticing that the function x, defined in Lemma 15.21 is given by 



X {y) = -nf3 - yn) 



one can check that is equivalent to (-ff^)- 

However, the function \ does not necessarily satisfy (Hi), which is too strong (as men- 
tioned in Remark [4| . For y > 0, 

a /rl 2nn„,n—l 

Consequently, \ is decreasing for y < 9 and increasing for y > 6. Taking y = f3~ 1 (5), we 
find that \ 1S decreasing on [0,/3 _1 (<5)] if and only if /3q < 25. In this case (Hi) is fulfilled. 
If (3o > 25, then \ is decreasing on the interval [0,9] and increasing on [9, P~ 1 (5)], yet r* is 
uniquely defined. 
Note that 

A(r)= 26 and B{r) = ^Poe^ - n5[(2e^ - l)p - 6] 

V ; 2e- dT -l V ' (2e- ST -l)P 
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Figure 1: The functions Zq (left) and Z\ (right) are drawn on the interval [0,t*) for param- 
eters given by (|3~4")) and n = 12. One can see that Zo has exactly two roots, t\ ~ 4.52 and 
T2 ~ 8.36, and Z\ has no root. 



with 

2S 2 e^ 2S 2 e- 5T (p + n5) 

(T) = {2e-*r-iy and B = { 2e-*T i)5ft • 

Then, assuming that (|3"3"|) holds true, we define the Zfc functions, as in ([2"5f . for t s [0,t*) 

by 

arccos ( J + 2fc7r 

We choose the parameters according to [2U [Ml ED] : 

5 = 0.05 days" 1 , f3 a = 1.77 days" 1 , 0=1. (34) 

Notice that the value of 9 is in fact normalized and does not influence the stability of system 
©-((J]) since all coefficients actually do not depend on 9. The value of 9 only influences the 
shape of the oscillations and the values of the steady states. 

Using Maple to determine the roots of Z n , we first check that Zq (and consequently all 
Zk functions) is strictly negative on [0, r*) for n < 10. Hence, from Theorem l5.ll the positive 
steady state E* = (S*,N*) of ©-® is locally asymptotically stable for r £ [0,t). 

For n > 10, Pujo-Menjouet et al. [29l [30] noticed, for the model ©-© with the in- 
troduction rate (3 depending only upon the nonprolifcrating phase population N(t), that 
oscillations may be observed. 

We choose n — 12, in keeping with values in [29, 30 . Then, we find that 

r ~ 13.3 days and t* ~ 9.66 days. 

One can see on Figure Q] that Zq has two positive roots in this case, T\ ~ 4.52 days and 
T2 ~ 8.36 days, and that Z\ is strictly negative, so all Z^ functions, with k > 1 have no roots. 
Consequently, there exist two critical values, t\ and T2, for which a stability switch can occur 
at E*. 

For r < T\ , one can check that the populations are asymptotically stable on Figure [U In 
this case r = 3.5 days and the solutions of (J3j) — (J4j) oscillate transiently to the steady state. 
Numerical simulations of the solutions of ©^Q are carried out with dde23 [3T], a Matlab 
solver for delay differential equations. 
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Figure 2: For r = 3.5 days, and the other parameters given by ([34]) with n = 12, the 
solutions S(t) (dashed line) and N(t) (solid line) oscillate transiently to the steady state, 
which is asymptotically stable. Damped oscillations are observed. 




(a) (b) 



Figure 3: For r = 4.52 days, and the other parameters given by (|3"4"]) with n = 12, a 
Hopf bifurcation occurs and the steady state (S*,N*) of ([3])-(j4]) is unstable. The periodic 
solutions S(t) (dashed line) and N(t) (solid line) are represented in (a), and we can observe 
the solutions in the (S, iV)-plane in (b). Periods of the oscillations are about 15 days. 



When t = T\ , one can check that 

B(A 2 +A' + AA't c ) - B 3 - B 2 B't c - B'A ~ 0.053, 

so condition (f32|) holds, and a Hopf bifurcation occurs at (S*,N*), from Theorem 15. II This 
is illustrated on Figure[3] Periodic solutions with periods about 15 days are observed at the 
bifurcation, and the steady state E* becomes unstable. 

When r increases after the bifurcation, one can observe oscillating solutions with longer 
periods (in the order of 20 to 30 days), as it can be seen in Figure HI 

This phenomenon has already been observed by Pujo-Menjouet et al. [29l [30] . It can 
be related to diseases affecting blood cells, the so-called periodic hematological diseases |19j . 
which are characterized by oscillations of circulating blood cell counts with long periods 
compared to the cell cycle duration. Among the wide variety of periodic hematological 
diseases, we can cite chronic myelogenous leukemia [H HSl I2S] , a cancer of white blood cells 
with periods usually falling in the range of 70 to 80 days, and cyclical neutropenia [TO] [19] 
which is known to exhibit oscillations around 3 weeks of circulating neutrophils (white cells), 
as observed on Figured] 
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Figure 5: For r = 9 days, and the other parameters given by (|34p with n — 12, damped 
oscillations are observed and the steady state is stable. 



Eventually, one can note that when r passes through the second critical value ti , stability 
switches and the steady state (S*, N*) becomes stable again (see Figured]). 



7 Discussion 

We considered a nonlinear model of blood cell dynamics in which the nonlinearity depends 
upon the entire hematopoietic stem cell population, contrary to the common assumption 
used in previous works [6] [TO] EH E2] EH [30] dealing with blood cell production models. 
Then we were lead to the study of a new nonlinear system of two differential equations with 
delay (describing the cell cycle duration) modelling the hematopoietic stem cells dynamics. 

We obtained the existence of two steady states for this model: a trivial one and a positive 
delay-dependent steady state. Through sections 0] and [5] we performed the stability analysis 
of our model. We determined necessary and sufficient conditions for the global asymptotic 
stability of the trivial steady state of system ((3j)— (|4]) , which describes the population's dying 
out. Using an approach proposed by Beretta and Kuang [S], we analyzed a first degree 
exponential polynomial characteristic equation with delay-dependent coefficients in order to 
obtain the existence of a Hopf bifurcation for the positive steady state (see Theorem 15.11) , 
leading to the existence of periodic solutions. 

On the example presented in the previous section, we obtained long periods oscillations, 
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which can be related to some periodic hematological diseases (in particularly, to cyclical 
neutropenia [10j). This result is in keeping with previous analysis of blood cell dynamics 
models (as it can be found in [5TJ [5^1 [3D] ) . Periodic hematological diseases are particular 
diseases mostly originated from the hematopoietic stem cell compartment. The appearance of 
periodic solutions in our model with periods that can be related to the ones observed in some 
periodic hematological diseases stresses the interesting properties displayed by our model. 
Periods of oscillating solutions can for example be used to determine the length of cell cycles 
in hematopoietic stem cell populations that cannot be directly determined experimentally. 

Moreover, stability switches have been observed, due to the structure of the equations 
(nonlinear equations with delay-dependent coefficients). Such a behavior had been noted in 
previous works dealing with blood cell production models (see [29j [30] ) , but it had never 
been mathematically explained. 

We can note that our assumption that proliferating and nonproliferating cells die with 
the same rate may be too limitative, since Pujo-Menjouet et al. [HI OS] already noticed that 
the apoptotic rate (the proliferating phase mortality rate 7) plays an important role in the 
appearance of oscillating solutions. However, by assuming that the two populations die with 
different rates, we are lead to a second order exponential polynomial characteristic equation, 
and the calculations are more difficult than the ones carried out in the present work. We let 
it for further analysis. 
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